{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a8b892c8",
   "metadata": {},
   "source": [
    "Printed and electronic copies of *Modeling and Simulation in Python* are available from [No Starch Press](https://nostarch.com/modeling-and-simulation-python) and [Bookshop.org](https://bookshop.org/p/books/modeling-and-simulation-in-python-allen-b-downey/17836697?ean=9781718502161) and [Amazon](https://amzn.to/3y9UxNb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "funded-utilization",
   "metadata": {},
   "source": [
    "# The Empire State Building Strikes Back"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "imported-table",
   "metadata": {
    "tags": []
   },
   "source": [
    "*Modeling and Simulation in Python*\n",
    "\n",
    "Copyright 2021 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "electoral-turkey",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# install Pint if necessary\n",
    "\n",
    "try:\n",
    "    from pint import UnitRegistry\n",
    "except ImportError:\n",
    "    !pip install pint\n",
    "    \n",
    "# import units\n",
    "from pint import UnitRegistry\n",
    "units = UnitRegistry()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "formal-context",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# download modsim.py if necessary\n",
    "\n",
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://raw.githubusercontent.com/AllenDowney/' +\n",
    "         'ModSimPy/master/modsim.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "progressive-typing",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import functions from modsim\n",
    "\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "removable-zoning",
   "metadata": {
    "tags": []
   },
   "source": [
    "This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. \n",
    "Click here to access the notebooks: <https://allendowney.github.io/ModSimPy/>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "embedded-gentleman",
   "metadata": {},
   "source": [
    "So far the differential equations we've worked with have been *first order*, which means they involve only first derivatives. In this\n",
    "chapter, we turn our attention to *second order* differential equations, which can involve both first and second derivatives.\n",
    "\n",
    "We'll revisit the falling penny example from Chapter 1, and use `run_solve_ivp` to find the position and velocity of the penny as it falls, with and without air resistance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "isolated-louis",
   "metadata": {},
   "source": [
    "## Newton's Second Law of Motion\n",
    "\n",
    "First order differential equations (DEs) can be written \n",
    "\n",
    "$$\\frac{dy}{dx} = G(x, y)$$ \n",
    "\n",
    "where $G$ is some function of $x$ and $y$ (see <http://modsimpy.com/ode>). Second order DEs can be written \n",
    "\n",
    "$$\\frac{d^2y}{dx^2} = H(x, y, \\frac{dy}{dt})$$\n",
    "\n",
    "where $H$ is a function of $x$, $y$, and $dy/dx$.\n",
    "\n",
    "In this chapter, we will work with one of the most famous and useful\n",
    "second order DEs, Newton's second law of motion: \n",
    "\n",
    "$$F = m a$$ \n",
    "\n",
    "where $F$ is a force or the total of a set of forces, $m$ is the mass of a moving object, and $a$ is its acceleration."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "drawn-symphony",
   "metadata": {},
   "source": [
    "Newton's law might not look like a differential equation, until we\n",
    "realize that acceleration, $a$, is the second derivative of position,\n",
    "$y$, with respect to time, $t$. With the substitution\n",
    "\n",
    "$$a = \\frac{d^2y}{dt^2}$$ \n",
    "\n",
    "Newton's law can be written\n",
    "\n",
    "$$\\frac{d^2y}{dt^2} = F / m$$ \n",
    "\n",
    "And that's definitely a second order DE.\n",
    "In general, $F$ can be a function of time, position, and velocity."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "swiss-vietnam",
   "metadata": {},
   "source": [
    "Of course, this \"law\" is really a model in the sense that it is a\n",
    "simplification of the real world. Although it is often approximately\n",
    "true:\n",
    "\n",
    "-   It only applies if $m$ is constant. If mass depends on time,\n",
    "    position, or velocity, we have to use a more general form of\n",
    "    Newton's law (see <http://modsimpy.com/varmass>).\n",
    "\n",
    "-   It is not a good model for very small things, which are better\n",
    "    described by another model, quantum mechanics.\n",
    "\n",
    "-   And it is not a good model for things moving very fast, which are\n",
    "    better described by yet another model, relativistic mechanics.\n",
    "\n",
    "However, for medium-sized things with constant mass, moving at\n",
    "medium-sized speeds, Newton's model is extremely useful. If we can\n",
    "quantify the forces that act on such an object, we can predict how it\n",
    "will move."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "coordinate-three",
   "metadata": {},
   "source": [
    "## Dropping Pennies\n",
    "\n",
    "As a first example, let's get back to the penny falling from the Empire State Building, which we considered in Chapter 1. We will implement two models of this system: first without air resistance, then with.\n",
    "\n",
    "Given that the Empire State Building is 381 m high, and assuming that\n",
    "the penny is dropped from a standstill, the initial conditions are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "compatible-increase",
   "metadata": {},
   "outputs": [],
   "source": [
    "init = State(y=381, v=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "intellectual-radiation",
   "metadata": {},
   "source": [
    "where `y` is height above the sidewalk and `v` is velocity. \n",
    "\n",
    "I'll put the initial conditions in a `System` object, along with the magnitude of acceleration due to gravity, `g`, and the duration of the simulations, `t_end`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "reverse-authorization",
   "metadata": {},
   "outputs": [],
   "source": [
    "system = System(init=init, \n",
    "                g=9.8, \n",
    "                t_end=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "heavy-boards",
   "metadata": {},
   "source": [
    "Now we need a slope function, and here's where things get tricky. As we have seen, `run_solve_ivp` can solve systems of first order DEs, but Newton's law is a second order DE. However, if we recognize that\n",
    "\n",
    "1.  Velocity, $v$, is the derivative of position, $dy/dt$, and\n",
    "\n",
    "2.  Acceleration, $a$, is the derivative of velocity, $dv/dt$,\n",
    "\n",
    "we can rewrite Newton's law as a system of first order ODEs:\n",
    "\n",
    "$$\\frac{dy}{dt} = v$$ \n",
    "\n",
    "$$\\frac{dv}{dt} = a$$ \n",
    "\n",
    "And we can translate those\n",
    "equations into a slope function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "occupied-mercury",
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func(t, state, system):\n",
    "    y, v = state\n",
    "\n",
    "    dydt = v\n",
    "    dvdt = -system.g\n",
    "    \n",
    "    return dydt, dvdt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "opening-adolescent",
   "metadata": {},
   "source": [
    "As usual, the parameters are a time stamp, a `State` object, and a `System` object.\n",
    "\n",
    "The first line unpacks the state variables, `y` and `v`.\n",
    "\n",
    "The next two lines compute the derivatives of the state variables, `dydt` and `dvdt`.\n",
    "The derivative of position is velocity, and the derivative of velocity is acceleration.\n",
    "In this case, $a = -g$, which indicates that acceleration due to gravity is in the direction of decreasing $y$. \n",
    "\n",
    "`slope_func` returns a sequence containing the two derivatives.\n",
    "\n",
    "Before calling `run_solve_ivp`, it is a good idea to test the slope\n",
    "function with the initial conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "positive-feeling",
   "metadata": {},
   "outputs": [],
   "source": [
    "dydt, dvdt = slope_func(0, system.init, system)\n",
    "dydt, dvdt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "false-charlotte",
   "metadata": {},
   "source": [
    "The result is 0 m/s for velocity and -9.8 m/s$^2$ for acceleration.\n",
    "\n",
    "Now we call `run_solve_ivp` like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "lovely-management",
   "metadata": {},
   "outputs": [],
   "source": [
    "results, details = run_solve_ivp(system, slope_func)\n",
    "details.message"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ranging-lingerie",
   "metadata": {},
   "source": [
    "`results` is a `TimeFrame` with two columns: `y` contains the height of the penny; `v` contains its velocity.\n",
    "Here are the first few rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "assisted-swimming",
   "metadata": {},
   "outputs": [],
   "source": [
    "results.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "solved-chambers",
   "metadata": {},
   "source": [
    "We can plot the results like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "authorized-barrier",
   "metadata": {},
   "outputs": [],
   "source": [
    "results.y.plot()\n",
    "\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Position (m)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "differential-airfare",
   "metadata": {},
   "source": [
    "Since acceleration is constant, velocity increases linearly and position decreases quadratically; as a result, the height curve is a parabola.\n",
    "\n",
    "The last value of `results.y` is negative, which means we ran the simulation too long. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "protected-fiber",
   "metadata": {},
   "outputs": [],
   "source": [
    "results.iloc[-1].y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "metallic-tamil",
   "metadata": {},
   "source": [
    "One way to solve this problem is to use the results to\n",
    "estimate the time when the penny hits the sidewalk.\n",
    "\n",
    "The ModSim library provides `crossings`, which takes a `TimeSeries` and a value, and returns a sequence of times when the series passes through the value. We can find the time when the height of the penny is `0` like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "japanese-clear",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_crossings = crossings(results.y, 0)\n",
    "t_crossings"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "demonstrated-emission",
   "metadata": {},
   "source": [
    "The result is an array with a single value, 8.818 s. Now, we could run\n",
    "the simulation again with `t_end = 8.818`, but there's a better way."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "blind-dominant",
   "metadata": {},
   "source": [
    "## Events\n",
    "\n",
    "As an option, `run_solve_ivp` can take an *event function*, which\n",
    "detects an \"event\", like the penny hitting the sidewalk, and ends the\n",
    "simulation.\n",
    "\n",
    "Event functions take the same parameters as slope functions, `t`, `state`, and `system`. They should return a value that passes through `0` when the event occurs. Here's an event function that detects the penny hitting the sidewalk:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "comfortable-simple",
   "metadata": {},
   "outputs": [],
   "source": [
    "def event_func(t, state, system):\n",
    "    y, v = state\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "closing-vehicle",
   "metadata": {},
   "source": [
    "The return value is the height of the penny, `y`, which passes through\n",
    "`0` when the penny hits the sidewalk.\n",
    "\n",
    "We pass the event function to `run_solve_ivp` like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "exotic-shareware",
   "metadata": {},
   "outputs": [],
   "source": [
    "results, details = run_solve_ivp(system, slope_func,\n",
    "                                 events=event_func)\n",
    "details.message"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "recreational-blair",
   "metadata": {},
   "source": [
    "Then we can get the flight time like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "appropriate-roberts",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_end = results.index[-1]\n",
    "t_end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "pediatric-portal",
   "metadata": {},
   "source": [
    "And the final velocity like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "orange-retro",
   "metadata": {},
   "outputs": [],
   "source": [
    "y, v = results.iloc[-1]\n",
    "y, v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cleared-jamaica",
   "metadata": {},
   "source": [
    "If there were no air resistance, the penny would hit the sidewalk (or someone's head) at about 86 m/s. So it's a good thing there is air resistance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "induced-albert",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "In this chapter, we wrote Newton's second law, which is a second order DE, as a system of first order DEs.\n",
    "Then we used `run_solve_ivp` to simulate a penny dropping from the Empire State Building in the absence of air resistance.\n",
    "And we used an event function to stop the simulation when the penny reaches the sidewalk.\n",
    "\n",
    "In the next chapter we'll add air resistance to the model.\n",
    "But first you might want to work on this exercise."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "operational-bhutan",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. \n",
    "You can access the notebooks at <https://allendowney.github.io/ModSimPy/>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "straight-johns",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Here's a question from the web site *Ask an Astronomer* (see http://curious.astro.cornell.edu/about-us/39-our-solar-system/the-earth/other-catastrophes/57-how-long-would-it-take-the-earth-to-fall-into-the-sun-intermediate):\n",
    "\n",
    "> \"If the Earth suddenly stopped orbiting the Sun, I know eventually it would be pulled in by the Sun's gravity and hit it. How long would it take the Earth to hit the Sun? I imagine it would go slowly at first and then pick up speed.\"\n",
    "\n",
    "Use `run_solve_ivp` to answer this question.\n",
    "\n",
    "Here are some suggestions about how to proceed:\n",
    "\n",
    "1.  Look up the Law of Universal Gravitation and any constants you need.  I suggest you work entirely in SI units: meters, kilograms, and Newtons.\n",
    "\n",
    "2.  When the distance between the Earth and the Sun gets small, this system behaves badly, so you should use an event function to stop when the surface of Earth reaches the surface of the Sun.\n",
    "\n",
    "3. Express your answer in days, and plot the results as millions of kilometers versus days.\n",
    "\n",
    "If you read the reply by Dave Rothstein, you will see other ways to solve the problem, and a good discussion of the modeling decisions behind them.\n",
    "\n",
    "You might also be interested to know that it's not that easy to get to the Sun; see https://www.theatlantic.com/science/archive/2018/08/parker-solar-probe-launch-nasa/567197/."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "forbidden-distributor",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "former-taxation",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "oriental-riverside",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "radio-reproduction",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "heavy-cologne",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "little-electric",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "continental-details",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "suitable-traveler",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "upper-victory",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "transparent-treat",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "brutal-woman",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "gentle-burst",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "comfortable-galaxy",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "satisfactory-latitude",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "significant-rebound",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "listed-shelter",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
